In this section, we develop a set of hierarchical Bayesian models to
analyze and explain the factors influencing worker productivity in the
garment sector. Given the complexity and nested nature of the data,
where individual observations are clustered within production teams, we
adopt a modeling framework that allows for both global effects and
team-specific variability.
The Bayesian paradigm offers several advantages in this context: it
enables the incorporation of prior knowledge, provides full posterior
distributions for uncertainty quantification and handles hierarchical
structures in a principled and coherent way. We begin with a relatively
simple centered hierarchical model, that already includes non-linear
effects through spline bases, and then progressively increase model
complexity by testing different interaction structures and modifying the
hierarchical specification. We explore both centered and non-centered
parameterizations for the random team effects, which allows us to assess
how reparameterization impacts convergence, identifiability and model
fit.
Through a combination of convergence diagnostics, posterior inference
and model comparison criteria such as the Deviance Information Criterion
(DIC) and posterior predictive checks, we systematically evaluate the
performance of each model. This stepwise approach highlights how each
component contributes to the model’s ability to capture the key drivers
of productivity and ultimately guides us toward a flexible and
interpretable specification that best fits the observed data.
Model 1: Centered
Hierarchical Additive Model
This first Bayesian model serves as the baseline and
was constructed by integrating the insights gained from both EDA and the
results of the GAMMs. The model is designed to combine linear effects,
smooth nonlinear effects and a hierarchical structure that captures
variability across teams.
The response variable \(y_i\), which
represents actual_productivity_nqt, is assumed to follow a
Student-t distribution:
\[
y_i \sim \text{Student-t}(\mu_i, \tau^{-1}, \nu)
\]
We use this distribution instead of a normal distribution because it
is more robust to extreme values. The degrees of freedom \(\nu\), which control the heaviness of the
tails, are estimated from the data using an exponential
prior:
\[
\nu \sim \text{Exponential}(1)
\]
This choice allows the model to flexibly adapt to the presence of
outliers: smaller values of \(\nu\)
imply heavier tails and greater robustness, while larger values
approximate a normal distribution. The prior favors moderate to low
values, letting the data inform the appropriate tail behavior.
The mean of the response \(y_i\) is
modeled as the sum of several components as follows
\[
\mu_i = \beta_0 + \mathbf{X}_i^\top \boldsymbol{\beta}_X +
\mathbf{Z}_{1,i}^\top \boldsymbol{\beta}_1 + \mathbf{Z}_{2,i}^\top
\boldsymbol{\beta}_2 + (u_{j[i]} - \bar{u}),
\]
where:
- \(\beta_0\) is the overall
intercept;
- \(\mathbf{X}_i\) contains the
linear predictors
targeted_productivity,
incentive and idle_time_flag;
- \(\mathbf{Z}_{1,i}\) and \(\mathbf{Z}_{2,i}\) contain spline basis
functions for
over_time and smv. These
variables showed nonlinear patterns in EDA, which we represent using
cubic regression splines;
- \(u_{j[i]}\) is a team-specific
intercept that accounts for unobserved differences between teams. We
subtract the mean \(\bar{u}\) to ensure
that the group effects are centered and do not absorb the global
intercept.
We chose to use orthogonal spline bases for the
nonlinear components, because it helps reduce multicollinearity between
the basis functions and improves the convergence of the Markov Chain
Monte Carlo (MCMC) algorithm. Orthogonal bases ensure that each
smooth term contributes uniquely to the model, making parameter
estimation more stable and efficient.
All regression coefficients are assigned standard normal
priors:
\[
\beta \sim \mathcal{N}(0, 1)
\]
These priors are weakly informative: they provide regularization,
which helps prevent overfitting, but they are not overly
restrictive.
For the standard deviation of the residuals, we use:
\[
\sigma \sim \text{Uniform}(0, 100)
\]
The team-level effects \(u_j\) are
modeled with a normal distribution:
\[
u_j \sim \mathcal{N}(0, \tau_{\text{team}}^{-1})
\]
The prior for the team-level standard deviation \(\text{sd_team}\) is a
half-Student-t distribution:
\[
\text{sd_team} \sim \text{Student-}t^+(0, 1, 1)
\]
This prior is commonly used in hierarchical models because it allows
for flexible group-level variation while providing some shrinkage, which
improves stability especially when the number of groups is small.
In summary, this model provides a flexible and
interpretable structure that balances
robustness, generalization and
computational stability. It captures the most relevant
features of the data and serves as a solid starting point for more
complex or specialized models.
Model Definition
and Setup
We now provide the complete definition of the Bayesian model and the
corresponding data structures used for inference. The model is specified
using JAGS syntax and reflects the formulation previously discussed in
the theoretical description. We also define the dataset structure passed
to JAGS, ensuring that the design matrices for the fixed and smooth
components are constructed using orthogonalized spline bases. This
improves identifiability and MCMC efficiency by reducing collinearity
between predictors. Initial values for the Markov chains are defined in
a way that supports model convergence, with small random perturbations
around the empirical mean for each parameter. In particular, the
team-level effects are initialized to have mean zero to preserve
identifiability of the overall intercept.
set.seed(123)
model_string <- "
model {
for (i in 1:N) {
mu[i] <- beta0 +
inprod(betaX[], X[i,]) +
inprod(beta1[], Z1[i,]) +
inprod(beta2[], Z2[i,]) +
(u[team[i]] - u_mean)
y[i] ~ dt(mu[i], tau, nu)
y_rep[i] ~ dt(mu[i], tau, nu)
}
sigma ~ dunif(0, 100)
tau <- pow(sigma, -2)
nu ~ dexp(1)
beta0 ~ dnorm(0, 1)
for (j in 1:KX) {
betaX[j] ~ dnorm(0, 1)
}
for (j in 1:K1) {
beta1[j] ~ dnorm(0, 1)
}
for (j in 1:K2) {
beta2[j] ~ dnorm(0, 1)
}
for (j in 1:J_team) {
u[j] ~ dnorm(0, tau_team)
}
sd_team ~ dt(0, 1, 1) T(0,)
tau_team <- pow(sd_team, -2)
u_mean <- mean(u[])
}
"
writeLines(model_string, con = "model.jags")
X <- model.matrix(~ 0 + targeted_productivity + incentive + idle_time_flag, data = data)
Z_overtime <- smoothCon(s(over_time, bs = "cr", k = 3), data = data, absorb.cons = TRUE)[[1]]$X
Z_smv <- smoothCon(s(smv, bs = "cr", k = 3), data = data, absorb.cons = TRUE)[[1]]$X
jags_data <- list(
y = as.vector(data$actual_productivity_nqt),
X = X,
Z1 = Z_overtime,
Z2 = Z_smv,
team = as.numeric(as.factor(data$team)),
N = nrow(data),
J_team = length(unique(data$team)),
KX = ncol(X),
K1 = ncol(Z_overtime),
K2 = ncol(Z_smv)
)
inits <- function() {
beta0_init <- mean(jags_data$y)
u_init <- rnorm(jags_data$J_team, 0, 0.1)
u_centered <- u_init - mean(u_init)
list(
beta0 = beta0_init,
betaX = rnorm(jags_data$KX, 0, 0.1),
beta1 = rnorm(jags_data$K1, 0, 0.1),
beta2 = rnorm(jags_data$K2, 0, 0.1),
u = u_centered,
sigma = runif(1, 0.1, 1),
sd_team = runif(1, 0.1, 1),
nu = runif(1, 2, 10)
)
}
Parameter Selection
and Model Sampling
We define the parameters to monitor during the MCMC sampling.
Diagnostic parameters include all regression coefficients, the random
effect standard deviation and the group-level effects. In addition, we
sample replicated responses y_rep for posterior predictive
checking and request the calculation of the Deviance Information
Criterion (DIC) for model evaluation.
params_diag <- c("betaX", "beta0", "beta1", "beta2", "u", "sd_team")
params_ppc <- c("y_rep")
jags_model <- jags(
data = jags_data,
inits = inits,
parameters.to.save = c(params_diag, params_ppc, "deviance"),
model.file = "model.jags",
n.chains = 3,
n.iter = 10000,
n.burnin = 2000,
n.thin = 1,
DIC = TRUE
)
samples_diag <- as.mcmc(jags_model)
To ensure reproducibility and speed up report rendering, we load the
precomputed posterior samples and model fit object from disk.
samples_diag <- readRDS("samples_diag_simple.rds")
jags_model <- readRDS("jags_model_simple.rds")
Posterior
Predictive Check (PPC)
We now evaluate how well the model replicates the observed data by
comparing the distribution of observed outcomes to the distribution of
simulated responses drawn from the posterior predictive
distribution.
The posterior predictive check reveals a generally good agreement
between the observed and predicted distributions of normalized actual
productivity. The central bulk of the distribution, particularly around
the mode near zero, is well captured by the model, suggesting that the
main trends in the data are adequately represented.
However, we observe two minor discrepancies:
- The predicted distribution appears slightly
smoother and more peaked near the central mode,
potentially indicating some degree of over-regularization.
- The left tail (around -2) is
underestimated by the model, suggesting that some
extreme low-productivity values are not fully captured.
- The right tail is slightly underpredicted, though
this deviation is relatively minor.
Overall, the model provides a satisfactory fit, with the posterior
predictive distribution broadly capturing the shape and spread of the
observed data. These results support the adequacy of the model structure
and prior specification, though minor tail mismatches could potentially
be addressed in future refinements.
Autocorrelation
Diagnostics
The autocorrelation plots displayed below provide insight into the
sampling efficiency of the MCMC algorithm for each parameter. For
readability, we report only the first 30 lags, which are generally
sufficient to assess short-term dependence in the chains.
Across all monitored parameters, autocorrelation declines rapidly
toward zero. This pattern indicates good mixing behavior and low
correlation between consecutive samples. Overall, the MCMC chains
exhibit desirable autocorrelation properties, which supports the
reliability of posterior summaries and the adequacy of the chosen priors
and model specification.
Traceplots of
Parameters and Posterior Density Overlays
We inspect traceplots of the MCMC chains to ensure proper mixing and
stationarity. Well-behaved chains should explore the posterior
distribution thoroughly without visible drift or trends.
In this case, all traceplots display stable horizontal
bands with no visible trends or drifts, which is a strong
indication of convergence. The chains for each
parameter show rapid oscillation around a constant
mean. Overall, the traceplots provide reassuring evidence that
the chains have reached stationarity and are adequately sampling from
the posterior distributions.
To further support our convergence assessment, we visualize the
posterior distributions of each parameter across the MCMC chains. The
high degree of overlap among the density curves confirms the consistency
of the chains and provides additional evidence of proper convergence in
this model.
Inference and
Parameter Interpretation
The table below presents summary statistics for all estimated
parameters, including posterior means, standard deviations, credible
intervals and diagnostic metrics such as R-hat and effective sample
sizes.
Posterior Summary
|
variable
|
mean
|
median
|
sd
|
mad
|
mcse_mean
|
mcse_sd
|
rhat
|
ess_bulk
|
ess_tail
|
2.5%
|
25%
|
50%
|
97.5%
|
|
beta0
|
-0.470
|
-0.470
|
0.020
|
0.021
|
0.001
|
0.000
|
1.007
|
1627.161
|
3717.687
|
-0.510
|
-0.484
|
-0.470
|
-0.430
|
|
beta1[1]
|
-0.129
|
-0.128
|
0.049
|
0.049
|
0.001
|
0.000
|
1.000
|
7610.198
|
11778.713
|
-0.226
|
-0.161
|
-0.128
|
-0.034
|
|
beta1[2]
|
-0.241
|
-0.239
|
0.070
|
0.070
|
0.001
|
0.001
|
1.001
|
5510.813
|
8032.367
|
-0.381
|
-0.287
|
-0.239
|
-0.107
|
|
beta2[1]
|
0.061
|
0.061
|
0.049
|
0.048
|
0.001
|
0.000
|
1.000
|
2316.974
|
4733.202
|
-0.036
|
0.029
|
0.061
|
0.155
|
|
beta2[2]
|
-0.981
|
-0.979
|
0.116
|
0.115
|
0.001
|
0.001
|
1.001
|
6904.462
|
12757.937
|
-1.210
|
-1.058
|
-0.979
|
-0.756
|
|
betaX[1]
|
0.248
|
0.248
|
0.014
|
0.013
|
0.000
|
0.000
|
1.003
|
2933.938
|
6799.576
|
0.221
|
0.239
|
0.248
|
0.275
|
|
betaX[2]
|
0.599
|
0.599
|
0.023
|
0.023
|
0.001
|
0.000
|
1.007
|
1530.500
|
3263.046
|
0.553
|
0.583
|
0.599
|
0.643
|
|
betaX[3]
|
-0.372
|
-0.366
|
0.157
|
0.155
|
0.001
|
0.001
|
1.000
|
12557.445
|
12169.283
|
-0.702
|
-0.473
|
-0.366
|
-0.080
|
|
sd_team
|
0.120
|
0.113
|
0.038
|
0.032
|
0.001
|
0.002
|
1.001
|
3122.536
|
5665.926
|
0.064
|
0.094
|
0.113
|
0.210
|
|
u[1]
|
0.094
|
0.092
|
0.053
|
0.052
|
0.001
|
0.001
|
1.001
|
2460.838
|
3889.824
|
-0.006
|
0.058
|
0.092
|
0.202
|
|
u[10]
|
-0.062
|
-0.061
|
0.049
|
0.047
|
0.001
|
0.001
|
1.001
|
2209.300
|
3611.534
|
-0.160
|
-0.093
|
-0.061
|
0.034
|
|
u[11]
|
-0.136
|
-0.134
|
0.061
|
0.061
|
0.001
|
0.001
|
1.001
|
1765.849
|
2872.654
|
-0.263
|
-0.176
|
-0.134
|
-0.023
|
|
u[12]
|
-0.117
|
-0.115
|
0.056
|
0.054
|
0.001
|
0.001
|
1.001
|
1732.252
|
2358.106
|
-0.236
|
-0.152
|
-0.115
|
-0.014
|
|
u[2]
|
0.155
|
0.153
|
0.059
|
0.058
|
0.001
|
0.001
|
1.001
|
3075.093
|
4917.095
|
0.047
|
0.115
|
0.153
|
0.277
|
|
u[3]
|
-0.052
|
-0.053
|
0.047
|
0.046
|
0.001
|
0.001
|
1.001
|
2101.652
|
3334.018
|
-0.144
|
-0.083
|
-0.053
|
0.042
|
|
u[4]
|
0.056
|
0.055
|
0.049
|
0.048
|
0.001
|
0.001
|
1.001
|
2227.494
|
4094.144
|
-0.039
|
0.023
|
0.055
|
0.154
|
|
u[5]
|
0.000
|
0.000
|
0.051
|
0.050
|
0.001
|
0.001
|
1.001
|
2422.792
|
4000.558
|
-0.099
|
-0.033
|
0.000
|
0.104
|
|
u[6]
|
-0.057
|
-0.055
|
0.051
|
0.051
|
0.001
|
0.001
|
1.002
|
1695.403
|
2725.233
|
-0.163
|
-0.090
|
-0.055
|
0.039
|
|
u[7]
|
0.106
|
0.104
|
0.053
|
0.052
|
0.001
|
0.001
|
1.001
|
2629.754
|
4557.211
|
0.008
|
0.070
|
0.104
|
0.213
|
|
u[8]
|
0.100
|
0.098
|
0.054
|
0.053
|
0.001
|
0.001
|
1.001
|
2640.975
|
4615.473
|
-0.004
|
0.063
|
0.098
|
0.211
|
|
u[9]
|
-0.088
|
-0.088
|
0.048
|
0.047
|
0.001
|
0.001
|
1.001
|
2079.430
|
3466.754
|
-0.187
|
-0.119
|
-0.088
|
0.005
|
All parameters exhibit excellent convergence diagnostics: R-hat
values are uniformly close to 1.000 and effective sample sizes largely
exceed 1000, ensuring reliable and stable posterior estimates.
The intercept (beta0) is significantly negative (mean =
–0.47, 95% CI: [–0.51, –0.43]) and, due to the centering of the
team-level effects, represents the average baseline productivity across
all teams after normalization. Among the linear predictors,
betaX[2] (incentive) shows a strong positive effect on
productivity (mean = 0.599, CI: [0.553, 0.643]) and
betaX[1] (likely targeted productivity) also exhibits a
positive association (mean = 0.248, CI: [0.221, 0.275]). In contrast,
betaX[3] (idle time flag) has a significant negative effect
(mean = –0.372, CI: [–0.702, –0.080]), indicating that flagged idle time
is associated with reduced productivity.
Spline coefficients for over_time (beta1)
and smv (beta2) capture flexible, non-linear
effects. These coefficients are not directly interpretable individually,
but their magnitude and joint significance indicate
how much and where the estimated function deviates from
linearity. While individual spline terms are not directly interpretable,
some coefficients (notably beta2[2] (mean ≈ –0.981),
beta1[1] and beta1[2]) have credible intervals
excluding zero, supporting the presence of non-linear relationships with
productivity.
The estimated standard deviation of the team-level effects
(sd_team) is around 0.12 and credibly positive, revealing
moderate variability between teams. In this model, team-level intercepts
u[j] have been centered, meaning each team effect is
expressed as a deviation from the mean team productivity. This
reparameterization improves model identifiability and allows the global
intercept to absorb the average effect. Most u[j] values
have intervals including zero, suggesting that many teams do not differ
substantially from the overall mean. However, u[2],
u[7], u[11] and u[12] have
credible intervals that exclude zero, indicating that these teams
exhibit statistically distinct behavior relative to the average.
The parameters with credible intervals that exclude zero and are
therefore considered statistically relevant include beta0,
betaX[1], betaX[2], betaX[3],
beta1[1], beta1[2], beta2[2],
sd_team and the team effects u[2],
u[7], u[11], u[12]. These
contribute most substantially to explaining variation in productivity.
For the remaining parameters, credible intervals contain zero,
reflecting greater posterior uncertainty about their individual impact,
though they may still contribute to capturing complex patterns or
improving regularization.
Model 2: Non-Centered
Hierarchical Model with Spline Interaction Incentive × WIP
Guided by our exploratory data analysis, we found strong evidence
that the variable wip may play an important nonlinear role
in determining productivity, especially through its interaction with
incentive and possible relation with smv.
While incentive was previously modeled linearly and
appeared statistically significant, we now test for its nonlinear
effects and interactions. We therefore propose a model that includes
both linear and smooth terms and a flexible interaction between
incentive and wip.
The outcome variable \(y_i\) is once
again modeled using a Student-t likelihood, improving
robustness against outliers by allowing heavier tails. Specifically:
\[
y_i \sim \text{Student-t}(\mu_i, \sigma, \nu)
\]
The linear predictor \(\mu_i\) is
composed as: \[
\mu_i = X_i^\top \beta_X + Z_{1,i}^\top \beta_1 + Z_{2,i}^\top \beta_2 +
Z_{3,i}^\top \beta_3 + u_{\text{j}[i]}
\] where:
- \(\beta_X\) are linear coefficients
for
targeted_productivity and idle_time_flag,
including also the intercept;
- \(\beta_1\), \(\beta_2\) and \(\beta_3\) are spline coefficients for
over_time, smv and the nonlinear interaction
between incentive and wip, respectively;
- \(u_{\text{team}[i]}\) is the
team-level random effect.
Smooth terms are constructed using orthogonalized
splines to ensure identifiability and better sampling
performance.
The group-level effects are modeled hierarchically:
- \(u_j = u_{\text{raw}, j} \cdot
\text{sd}_{\text{team}}\), with \(u_{\text{raw}, j} \sim \mathcal{N}(0,
1)\);
- \(\text{sd}_{\text{team}} \sim
\text{Uniform}(0, 100)\), allowing for flexibility in group-level
variability.
We deliberately do not center the group-level
effects in this model. As a result, the global intercept
betaX[1] represents the overall baseline productivity
across all teams, rather than being adjusted to reflect a mean-zero
constraint on team effects. This modeling choice allows each
u[j] to be interpreted as the absolute deviation of team
j from the global baseline, which is especially valuable
when the goal is to compare team performance directly.
The team-level effects are modeled using a non-centered
parameterization, where each u[j] is expressed as
the product of a standard normal variable and a shared standard
deviation parameter sd_team. This implies the prior \(u_j \sim \mathcal{N}(0, \text{sd_team}^2)\)
with sd_team following a weakly informative uniform prior,
allowing the data to inform how heterogeneous the teams are. The
non-centered approach enhances sampling efficiency and numerical
stability, particularly when the variance across teams is low or only
weakly informed by the data. Overall, this parameterization makes the
estimation process more stable and the interpretation of team effects
more transparent. In the previous model, where team effects were
centered, most team-specific estimates hovered around zero with few
statistically distinguishable deviations. By removing this constraint,
we aim to explore whether more meaningful absolute differences emerge
when each team’s deviation from the global baseline is modeled
directly.
To maintain identifiability, the prior on the intercept \(\beta_X[1] \sim \mathcal{N}(0, 0.1)\) is
tighter, which regularizes its value and prevents it from absorbing
excessive variance. All other coefficients, including splines and
interaction terms, receive weakly informative priors \(\mathcal{N}(0, 1)\), promoting flexibility
while avoiding overfitting.
The residual standard deviation \(\sigma\) is assigned a broad prior \(\text{Uniform}(0, 100)\) and the degrees of
freedom \(\nu\) of the Student-t
distribution are given an exponential prior \(\nu \sim \text{Exponential}(1)\).
In summary, this model combines additive splines, interaction terms
and random effects within a Student-t framework to provide a robust,
interpretable and flexible modeling approach informed by both prior
knowledge and data-driven structure.
Model Definition
and Setup
We now present the full specification of the extended Bayesian model,
along with the corresponding data and initialization setup.
set.seed(123)
model_string <- "
model {
for (i in 1:N) {
mu[i] <- inprod(betaX[], X[i,]) +
inprod(beta1[], Z1[i,]) +
inprod(beta2[], Z2[i,]) +
inprod(beta3[], Z3[i,]) +
u[team[i]]
y[i] ~ dt(mu[i], tau, nu)
y_rep[i] ~ dt(mu[i], tau, nu)
}
sigma ~ dunif(0, 100)
tau <- pow(sigma, -2)
nu ~ dexp(1)
betaX[1] ~ dnorm(0, 0.1)
for (j in 2:KX) {
betaX[j] ~ dnorm(0, 1)
}
for (j in 1:K1) {
beta1[j] ~ dnorm(0, 1)
}
for (j in 1:K2) {
beta2[j] ~ dnorm(0, 1)
}
for (j in 1:K3) {
beta3[j] ~ dnorm(0, 1)
}
sd_team ~ dunif(0, 100)
for (j in 1:J_team) {
u_raw[j] ~ dnorm(0, 1)
u[j] <- u_raw[j] * sd_team
}
tau_team <- pow(sd_team, -2)
}
"
writeLines(model_string, con = "model2.jags")
X <- model.matrix(~ targeted_productivity + idle_time_flag, data = data)
Z_overtime <- smoothCon(s(over_time, bs = "cr", k = 3), data = data, absorb.cons = TRUE)[[1]]$X
Z_smv <- smoothCon(s(smv, bs = "cr", k = 3), data = data, absorb.cons = TRUE)[[1]]$X
Z_interaction <- smoothCon(te(incentive, wip, bs = "cr", k = 3), data = data, absorb.cons = TRUE)[[1]]$X
jags_data <- list(
y = as.vector(data$actual_productivity_nqt),
X = X,
Z1 = Z_overtime,
Z2 = Z_smv,
Z3 = Z_interaction,
team = as.numeric(as.factor(data$team)),
N = nrow(data),
J_team = length(unique(data$team)),
KX = ncol(X),
K1 = ncol(Z_overtime),
K2 = ncol(Z_smv),
K3 = ncol(Z_interaction)
)
inits <- function() {
list(
betaX = rnorm(jags_data$KX, 0, 0.1),
beta1 = rnorm(jags_data$K1, 0, 0.1),
beta2 = rnorm(jags_data$K2, 0, 0.1),
beta3 = rnorm(jags_data$K3, 0, 0.1),
u_raw = rnorm(jags_data$J_team, 0, 1),
sigma = runif(1, 0.1, 1),
sd_team = runif(1, 0.1, 1),
nu = runif(1, 2, 10)
)
}
Parameter Selection
and Model Sampling
We monitor the same key parameters as in the previous model: fixed
effects, spline coefficients, team-level effects and their standard
deviation. We also generate replicated responses (y_rep)
for posterior predictive checks and compute the Deviance Information
Criterion (DIC) to assess model fit.
params_diag <- c("betaX", "beta1", "beta2", "beta3", "u", "sd_team")
params_ppc <- c("y_rep")
jags_model <- jags(
data = jags_data,
inits = inits,
parameters.to.save = c(params_diag, params_ppc, "deviance"),
model.file = "model2.jags",
n.chains = 3,
n.iter = 10000,
n.burnin = 2000,
n.thin = 1,
DIC = TRUE
)
samples_diag <- as.mcmc(jags_model)
As in the previous model, we save both the fitted model and MCMC
diagnostics using saveRDS() for reproducibility and future
analysis.
samples_diag <- readRDS("samples_diag2_simple.rds")
jags_model <- readRDS("jags_model2_simple.rds")
Posterior
Predictive Check (PPC)
We assess model fit by comparing the observed data distribution with
that of the posterior predictive simulations.
The PPC for Model 2 shows a substantial improvement over Model 1. The
predicted distribution more closely aligns with the observed
data, particularly in the central peak and the right tail.
While some discrepancies remain in the left tail and around local modes,
the model now captures the overall shape of the data
much more faithfully.
This suggests that the introduction of the spline-based interaction
between incentive and wip significantly
improved the model’s ability to replicate the empirical distribution of
productivity scores. The better alignment between predicted and observed
densities confirms that Model 2 provides a more realistic generative
process compared to the baseline specification in Model 1.
Autocorrelation
Diagnostics
We now report the autocorrelation diagnostics for Model 2, focusing
on the first 30 lags to evaluate sampling efficiency and potential
issues with chain mixing.
Overall, the autocorrelation drops off quickly for most parameters,
indicating good mixing and effective sampling. The decay is particularly
steep for the fixed effects and the random effect standard deviation
sd_team, confirming that the non-centered parameterization
continues to improve convergence properties. No persistent
autocorrelation is observed, suggesting that thinning is not required
and that the number of effective samples is adequate across monitored
parameters.
Traceplots of
Parameters and Posterior Density Overlays
To further assess MCMC convergence and sampling efficiency, we
inspect both the traceplots and the posterior density overlays for the
monitored parameters.
The traceplots reveal good mixing behavior across all chains, with no
evident trends or drifts. The chains appear stationary and
well-overlapped for each parameter, suggesting that the sampling process
has sufficiently explored the posterior distribution. Even the
group-level effects u[j] show stable and consistent traces,
supporting the effectiveness of the non-centered parameterization
adopted in the model.
The posterior density overlays confirm that the marginal posterior
distributions from each chain align almost perfectly. The densities
exhibit strong overlap and smooth unimodal shapes, which further
indicates that all chains are sampling from the same distribution.
Inference and
Parameter Interpretation
To assess the stability and uncertainty of the estimated parameters,
we present a comprehensive summary of their posterior distributions
alongside graphical representations of their 95% credible intervals.
Posterior Summary
|
variable
|
mean
|
median
|
sd
|
mad
|
mcse_mean
|
mcse_sd
|
rhat
|
ess_bulk
|
ess_tail
|
2.5%
|
25%
|
50%
|
75%
|
97.5%
|
|
beta1[1]
|
-0.093
|
-0.092
|
0.035
|
0.034
|
0.000
|
0.000
|
1.001
|
6595.368
|
9999.419
|
-0.164
|
-0.116
|
-0.092
|
-0.069
|
-0.026
|
|
beta1[2]
|
-0.037
|
-0.037
|
0.049
|
0.048
|
0.001
|
0.000
|
1.000
|
6086.241
|
9267.051
|
-0.132
|
-0.070
|
-0.037
|
-0.005
|
0.059
|
|
beta2[1]
|
0.069
|
0.069
|
0.032
|
0.032
|
0.001
|
0.000
|
1.001
|
2847.908
|
5595.702
|
0.004
|
0.047
|
0.069
|
0.090
|
0.131
|
|
beta2[2]
|
-0.763
|
-0.757
|
0.162
|
0.149
|
0.002
|
0.002
|
1.000
|
4667.709
|
6538.316
|
-1.121
|
-0.858
|
-0.757
|
-0.656
|
-0.461
|
|
beta3[1]
|
-0.185
|
-0.183
|
0.062
|
0.060
|
0.001
|
0.001
|
1.002
|
2143.030
|
3474.968
|
-0.313
|
-0.224
|
-0.183
|
-0.143
|
-0.068
|
|
beta3[2]
|
-0.459
|
-0.439
|
0.157
|
0.131
|
0.002
|
0.003
|
1.001
|
4879.730
|
5313.212
|
-0.835
|
-0.535
|
-0.439
|
-0.356
|
-0.213
|
|
beta3[3]
|
-0.140
|
-0.141
|
0.047
|
0.047
|
0.001
|
0.000
|
1.000
|
4989.977
|
8729.665
|
-0.233
|
-0.172
|
-0.141
|
-0.109
|
-0.046
|
|
beta3[4]
|
0.366
|
0.364
|
0.104
|
0.099
|
0.002
|
0.001
|
1.002
|
2620.341
|
4690.300
|
0.168
|
0.298
|
0.364
|
0.431
|
0.580
|
|
beta3[5]
|
-0.074
|
-0.075
|
0.056
|
0.056
|
0.001
|
0.000
|
1.001
|
5362.729
|
9130.119
|
-0.183
|
-0.112
|
-0.075
|
-0.036
|
0.039
|
|
beta3[6]
|
-0.046
|
-0.032
|
0.484
|
0.478
|
0.008
|
0.004
|
1.000
|
4127.047
|
7788.270
|
-1.030
|
-0.362
|
-0.032
|
0.280
|
0.868
|
|
beta3[7]
|
1.791
|
1.790
|
0.062
|
0.061
|
0.001
|
0.001
|
1.001
|
3824.381
|
6858.734
|
1.671
|
1.749
|
1.790
|
1.832
|
1.914
|
|
beta3[8]
|
2.021
|
2.027
|
0.153
|
0.148
|
0.003
|
0.002
|
1.001
|
3031.710
|
4471.204
|
1.702
|
1.925
|
2.027
|
2.125
|
2.305
|
|
betaX[1]
|
0.022
|
0.022
|
0.019
|
0.017
|
0.001
|
0.000
|
1.002
|
1220.384
|
2196.566
|
-0.014
|
0.011
|
0.022
|
0.033
|
0.060
|
|
betaX[2]
|
0.320
|
0.320
|
0.011
|
0.011
|
0.000
|
0.000
|
1.001
|
3413.679
|
5563.827
|
0.298
|
0.313
|
0.320
|
0.328
|
0.341
|
|
betaX[3]
|
-0.258
|
-0.224
|
0.184
|
0.136
|
0.002
|
0.002
|
1.000
|
7918.712
|
8162.251
|
-0.754
|
-0.326
|
-0.224
|
-0.140
|
0.011
|
|
sd_team
|
0.054
|
0.051
|
0.019
|
0.017
|
0.000
|
0.000
|
1.002
|
1619.544
|
2350.143
|
0.025
|
0.041
|
0.051
|
0.065
|
0.099
|
|
u[1]
|
0.026
|
0.025
|
0.034
|
0.033
|
0.001
|
0.000
|
1.001
|
3483.075
|
6496.685
|
-0.038
|
0.004
|
0.025
|
0.048
|
0.096
|
|
u[10]
|
-0.050
|
-0.049
|
0.028
|
0.027
|
0.001
|
0.000
|
1.001
|
2871.370
|
4984.004
|
-0.108
|
-0.067
|
-0.049
|
-0.030
|
0.003
|
|
u[11]
|
-0.038
|
-0.036
|
0.034
|
0.033
|
0.001
|
0.000
|
1.001
|
2629.598
|
4520.949
|
-0.108
|
-0.059
|
-0.036
|
-0.014
|
0.024
|
|
u[12]
|
-0.018
|
-0.017
|
0.031
|
0.030
|
0.001
|
0.000
|
1.001
|
2310.847
|
3525.315
|
-0.085
|
-0.038
|
-0.017
|
0.003
|
0.039
|
|
u[2]
|
0.025
|
0.024
|
0.029
|
0.028
|
0.001
|
0.000
|
1.001
|
2653.529
|
5007.223
|
-0.030
|
0.005
|
0.024
|
0.044
|
0.086
|
|
u[3]
|
-0.028
|
-0.027
|
0.027
|
0.027
|
0.001
|
0.000
|
1.001
|
2501.063
|
4806.296
|
-0.084
|
-0.046
|
-0.027
|
-0.010
|
0.024
|
|
u[4]
|
0.066
|
0.064
|
0.034
|
0.034
|
0.001
|
0.000
|
1.000
|
3238.721
|
5695.006
|
0.004
|
0.042
|
0.064
|
0.088
|
0.137
|
|
u[5]
|
-0.010
|
-0.011
|
0.030
|
0.028
|
0.001
|
0.000
|
1.001
|
3399.344
|
5705.707
|
-0.069
|
-0.029
|
-0.011
|
0.009
|
0.049
|
|
u[6]
|
-0.010
|
-0.009
|
0.028
|
0.026
|
0.001
|
0.000
|
1.001
|
2374.906
|
3774.745
|
-0.067
|
-0.027
|
-0.009
|
0.008
|
0.041
|
|
u[7]
|
0.054
|
0.053
|
0.034
|
0.033
|
0.001
|
0.000
|
1.000
|
3124.081
|
5986.409
|
-0.006
|
0.031
|
0.053
|
0.076
|
0.127
|
|
u[8]
|
0.040
|
0.038
|
0.032
|
0.031
|
0.001
|
0.000
|
1.000
|
3023.314
|
5615.285
|
-0.019
|
0.018
|
0.038
|
0.060
|
0.109
|
|
u[9]
|
-0.056
|
-0.054
|
0.027
|
0.027
|
0.001
|
0.000
|
1.001
|
2486.602
|
4443.467
|
-0.114
|
-0.073
|
-0.054
|
-0.037
|
-0.006
|
All parameters display excellent convergence properties: R-hat values
are tightly concentrated around 1.000 and effective sample sizes (both
bulk and tail) are consistently high often exceeding 3000, indicating
well-mixed chains and reliable posterior estimates.
The intercept (betaX[1]) is small and centered around
zero (mean = 0.022, 95% CI: [–0.014, 0.060]), as expected in a model
where the random team-level effects are not constrained to sum to zero.
This reflects the fact that the global baseline is absorbed through the
intercept without centering.
Among the linear predictors, betaX[2] which encodes
targeted_productivity has a strong positive effect (mean =
0.320, 95% CI: [0.298, 0.341]), suggesting a stable influence over time.
In contrast, betaX[3] (idle_time_flag)
displays a broad credible interval containing zero (mean = –0.258, CI:
[–0.754, 0.011]), indicating weak evidence for its effect.
The spline coefficients for over_time
(beta1), smv (beta2) and the
interaction between wip and incentive
(beta3) capture complex non-linear relationships. As with
all spline terms, the individual coefficients are not directly
interpretable in isolation, but their patterns and magnitudes offer
insight. Several terms have posterior means far from zero with narrow
credible intervals. Notably, beta3[7] and
beta3[8] show strong positive effects (means ≈ 1.79 and
2.02, respectively), indicating a strong non-linear contribution of the
interaction term to productivity. Also, beta3[4] and
beta2[1] exhibit a moderate positive effect, with posterior
means of approximately 0.366 and beta2[1] respectively.
Likewise, beta2[2] has a large negative mean (–0.763, CI:
[–1.121, –0.461]), revealing an important dip in the functional
relationship between smv and productivity. This suggests
that while smv may have a slight positive effect in certain
regions, the overall influence is predominantly negative. Moreover,
beta1[1] shows a slight negative effect (posterior mean =
–0.093), suggesting that this predictor may have a minor negative
influence on the target variable. Finally, beta3[1],
beta3[2] and beta3[3] show negative effects,
with posterior means of –0.185, –0.459 and –0.140, respectively. This
suggests that the interaction term may exhibit heterogeneous behavior:
in some regions it appears positive, in others negative and in others
not significant.
Interestingly, the spline coefficient beta3[6] has a
very wide credible interval (mean = –0.046, 95% CI: [–1.030, 0.868]),
suggesting this component may be unidentifiable or weakly supported by
the data, possibly due to multicollinearity or limited variation in the
spline basis at that location.
The estimated standard deviation of the team effects
(sd_team) is moderate and tightly estimated (mean = 0.054,
95% CI: [0.025, 0.099]), confirming the presence of non-negligible
heterogeneity among teams.
The team-specific intercepts u[j] are not centered and
should be interpreted as absolute deviations from the global baseline.
Several team effects exhibit credible intervals that exclude zero,
indicating statistically meaningful deviations in productivity. In
particular, u[4] (mean = 0.066, CI: [0.004, 0.137]) and
u[7] (mean = 0.054, CI: [–0.006, 0.127]) suggest
above-average productivity, while u[9] (mean = –0.056, CI:
[–0.114, –0.006]) reflects a significant underperformance. These
differences persist even after controlling for all included
covariates.
Parameters with credible intervals that exclude zero and are
therefore deemed statistically relevant include:
- Linear and spline terms betaX[2], beta1[1],
beta2[1], beta2[2], beta3[1],
beta3[2], beta3[3], beta3[4],
beta3[7], beta3[8]
- Group-level standard deviation sd_team
- Team effects u[4] and u[9]
These effects contribute most substantially to explaining
productivity differences. The other parameters have higher posterior
uncertainty and should be interpreted with caution, although they may
still play important roles in shaping smooth response surfaces or
improving model flexibility and predictive performance.
Model 3: Non-Centered
Hierarchical Model with Non-Linear Main Effect Incentive and Spline
Interaction WIP × SMV
This third model builds on the structure and insights gained from
Models 1 and 2 and aims to further refine the representation of
nonlinear relationships by modifying the specification of the smooth
terms.
As in the previous models, we use a Student-t distribution for the
observed outcomes, which is robust to outliers and heavy tails: \[
y[i] \sim \text{Student-t}(\mu[i], \tau^{-1}, \nu)
\] The predictor for each observation \(i\) is defined as: \[
\mu_i = \mathbf{X}_i^\top \boldsymbol{\beta}_X + \mathbf{Z}_{1,i}^\top
\boldsymbol{\beta}_1 + \mathbf{Z}_{2,i}^\top \boldsymbol{\beta}_2 +
\mathbf{Z}_{3,i}^\top \boldsymbol{\beta}_3 + u_{j[i]}
\]
Here, the linear terms X include
targeted_productivity and idle_time_flag,
while the smooth components consist of a univariate spline
Z1 for over_time, a bivariate tensor-product
spline Z2 for the interaction between smv and
wip and a univariate spline Z3 for
incentive. This choice is motivated by our exploratory data
analysis, which revealed that the interaction between smv
and wip may exhibit complex nonlinear effects on
productivity. While Model 2 previously included a spline-based
interaction between incentive and wipwhich led
to a substantial improvement in both fit and DIC, we observed that
incentive alone also acts as a strong linear predictor.
Nonetheless, we chose to retain a smooth (nonlinear) effect for
incentive in the current specification for two main
reasons. First, the prior results from Model 2 clearly indicate that
incentive can exhibit nonlinear behavior when
interacting with other variables, suggesting that its marginal
effect may not be entirely captured by a purely linear term. Second,
using a spline allows the model to flexibly adapt to any
plateaus, thresholds or saturation effects that a
linear term would fail to represent. Importantly, including
incentive as a univariate smooth term enables us to
disentangle the monetary effect from physical
production parameters such as smv and wip,
while still maintaining sufficient flexibility to capture its
nonlinearity. This design balances interpretability with model
performance and avoids the potential risk of overfitting associated with
more complex interactions.
The prior for the intercept betaX[1] is set as \(\mathcal{N}(0, 0.01^{-1})\), corresponding
to a variance of 100. This is a relatively weak prior that favors
stability while allowing the intercept to absorb the global baseline
across all teams. It is particularly important in our non-centered
framework, where the random effects are not constrained to sum to zero.
As a result, the intercept must capture the average productivity level
across all observations. The remaining regression coefficients and
spline weights betaX[2], beta1,
beta2, beta3 are assigned standard normal
priors \(\mathcal{N}(0,1)\), which
provide mild regularization without being overly restrictive. We use
broad \(\text{Uniform}(0, 100)\) priors
for the residual scale parameter sigma and the team-level
standard deviation sd_team, reflecting minimal prior
information. The degrees of freedom nu of the
t-distribution follow an \(\text{Exponential}(1)\) prior, which places
more mass on small values and therefore allows for heavier tails,
promoting robustness to outliers.
Team effects are modeled using a non-centered parameterization: \[
u[j] = u_{\text{raw}}[j] \cdot sd_{\text{team}} \quad {\text{where }}
u_{\text{raw}}[j] \sim \mathcal{N}(0, 1),
\] which we previously found to improve interpretability.
Model Definition
and Setup
We now define the JAGS model, the list of monitored parameters and
the function used to initialize the chains. The initialization is
designed to promote convergence and stability during sampling, while
remaining weakly informative.
In particular, we fix the intercept betaX[1] to a
moderately positive value (0.25), which reflects the average
productivity level observed in the data after standardization and
provides a sensible starting point in our non-centered framework. The
remaining regression coefficients and spline weights are initialized
from a narrow normal distribution centered at zero, allowing for slight
variation while avoiding extreme starting values. The team-level random
effects u_raw are sampled from a standard normal
distribution, consistent with their prior. The scale parameters
sigma and sd_team are initialized from a
uniform distribution on (0.1, 1), ensuring positivity and moderate
variability. Finally, the degrees of freedom nu of the
t-distribution are initialized from a uniform range between 2 and 10,
which avoids extreme tail behavior at the beginning of sampling.
set.seed(123)
model_string <- "
model {
for (i in 1:N) {
mu[i] <- inprod(betaX[], X[i,]) +
inprod(beta1[], Z1[i,]) +
inprod(beta2[], Z2[i,]) +
inprod(beta3[], Z3[i,]) +
u[team[i]]
y[i] ~ dt(mu[i], tau, nu)
y_rep[i] ~ dt(mu[i], tau, nu)
}
sigma ~ dunif(0, 100)
tau <- pow(sigma, -2)
nu ~ dexp(1)
betaX[1] ~ dnorm(0, 0.01)
for (j in 2:KX) {
betaX[j] ~ dnorm(0, 1)
}
for (j in 1:K1) {
beta1[j] ~ dnorm(0, 1)
}
for (j in 1:K2) {
beta2[j] ~ dnorm(0, 1)
}
for (j in 1:K3) {
beta3[j] ~ dnorm(0, 1)
}
sd_team ~ dunif(0, 100)
for (j in 1:J_team) {
u_raw[j] ~ dnorm(0, 1)
u[j] <- u_raw[j] * sd_team
}
tau_team <- pow(sd_team, -2)
}
"
writeLines(model_string, con = "model3.jags")
X <- model.matrix(~ targeted_productivity + idle_time_flag, data = data)
Z_overtime <- smoothCon(s(over_time, bs = "cr", k = 3), data = data, absorb.cons = TRUE)[[1]]$X
Z_incentive <- smoothCon(s(incentive, bs = "cr", k = 3), data = data, absorb.cons = TRUE)[[1]]$X
Z_interaction <- smoothCon(te(smv, wip, bs = "cr", k = 3), data = data, absorb.cons = TRUE)[[1]]$X
jags_data <- list(
y = as.vector(data$actual_productivity_nqt),
X = X,
Z1 = Z_overtime,
Z2 = Z_interaction,
Z3 = Z_incentive,
team = as.numeric(as.factor(data$team)),
N = nrow(data),
J_team = length(unique(data$team)),
KX = ncol(X),
K1 = ncol(Z_overtime),
K2 = ncol(Z_interaction),
K3 = ncol(Z_incentive)
)
inits <- function() {
list(
betaX = c(0.25, rnorm(jags_data$KX - 1, 0, 0.1)),
beta1 = rnorm(jags_data$K1, 0, 0.1),
beta2 = rnorm(jags_data$K2, 0, 0.1),
beta3 = rnorm(jags_data$K3, 0, 0.1),
u_raw = rnorm(jags_data$J_team, 0, 1),
sigma = runif(1, 0.1, 1),
sd_team = runif(1, 0.1, 1),
nu = runif(1, 2, 10)
)
}
Parameter Selection
and Model Sampling
We monitor the same categories of parameters as before: regression
terms, spline coefficients and group-level effects. We also include
replicated outcomes (y_rep) for posterior predictive checks
and the DIC for model comparison.
params_diag <- c("betaX", "beta1", "beta2", "beta3", "u", "sd_team")
params_ppc <- c("y_rep")
jags_model <- jags(
data = jags_data,
inits = inits,
parameters.to.save = c(params_diag, params_ppc, "deviance"),
model.file = "model3.jags",
n.chains = 3,
n.iter = 10000,
n.burnin = 2000,
n.thin = 1,
DIC = TRUE
)
As before, we save both the model object and the diagnostics with
saveRDS() for future access.
samples_diag <- readRDS("samples_diag3_simple.rds")
jags_model <- readRDS("jags_model3_simple.rds")
Posterior
Predictive Check (PPC)
We now perform a posterior predictive check to visually compare the
observed and predicted densities of the outcome variable.
The plot below shows a close alignment between the
predicted and observed densities, particularly around the main mode of
the distribution. Compared to the earlier models, this fit is visibly
improved: the underestimation in the left tail observed in Model 1 has
been largely corrected and the excessive smoothing of the central peak
seen in Model 2 has been mitigated.
The predicted distribution still slightly overshoots the observed one
in the region just below the mode, but the fit is clearly more faithful
overall. The improved match across the entire range, including in the
tails, confirms that this model captures both the
central tendency and dispersion more
effectively than previous specifications.
Autocorrelation
Diagnostics
We now examine the autocorrelation of the sampled parameters to
assess the efficiency of the MCMC chains.
The autocorrelation decays rapidly for all monitored parameters,
suggesting efficient sampling. Compared to previous models, the overall
structure indicates similar or slightly better mixing, particularly for
the spline coefficients and team-level effects.
Traceplots of
Parameters and Posterior Density Overlays
We assess the quality of MCMC convergence by examining both trace
plots and posterior density overlays for all monitored parameters.
The histograms reveal several key insights:
The boxplots display the distribution of standardized numerical
variables across the five quarters in the dataset. Overall, most
variables do not exhibit strong seasonal trends or significant shifts in
central tendency over time. However, a few subtle patterns can be
observed. For instance, 
